% In this program, we solve the FOC as written, allowing MatLab to
% determine the regime sequence

clear all
close all

global parms gamma

% The vector parms will contain constants and parameters that are not adjusted
% during the iterations. The parameters that are adjusted are:
% Alpha0 Alpha1 Alpha2 Alpha3 Gamma T2 N_T1 lambda_T1 k_T1 j_T1 k_int


y0 = 1.0;     % output of goods y0 - choose units so this is 1.0,

Sbar = 2126.0527;    % Total feasible resources
Res_0 = 1./0.065100642;       % Initial estimate of proved reserves
rho_0 = Sbar - Res_0;    % Initial ratio of Alpha2 to Alpha3
g_0 = 0.11068716;      % Initial value of per unit mining cost
n0 = 0.00832913;      % Initial level of mining investment

% gN_0 = -0.081639946;    % Partial derivative of g with respect to N at t=0
% Alpha3 = -gS_0.*Rho_0./gN_0;

alpha3 = 15; %-gS_0.*Rho_0./gN_0;
alpha2 = rho_0.*alpha3;

beta  =  0.05;
delta = 0.04;

c0 = 0.6619974; % Initial value of consumption for calculate lambda_0
k0 = 3.6071282734; % Initial value of K for the differential equation

A =  y0./k0;

% Marginal cost of backstop energy p = (Gamma1+H)^(-alpha)
% Alpha is the slope of the learning curve

alpha = 0.25;

% Given alpha, Gamma1 determines the initial cost of renewable energy. Here we set it
% to 4 times the initial cost of fossil fuel.

Gamma1 = (4*g_0).^(-1./alpha);

% After some time t, marginal cost will decline to Gamma2 and remain there.
% We assume this ultimate minimum marginal cost of renewables
% is 20% of the initial cost p when H = 0:

Gamma2 =0.8*g_0; 

% At Tbar, the marginal cost of backstop technology becomes Gamma2. Hence
% we have (Gamma1+H)^(-alpha)=Gamma2

H_Tbar = Gamma2^(-1/alpha)-Gamma1;

% Abar = beta-A*(1-Gamma2)+ delta
% Abar is a constant in solution of terminal regime differential equation
% used in the solutions for lambda and k.
% Note: We should have Abar < 0

Abar = beta+delta-A*(1-Gamma2); 

% psi is the effect of learning by doing on H
% psi has to be between 0 and 1
% A smaller value of psi will allow a larger role for explicit investment
% in renewable technology as opposed to learning by doing

psi = 0.33;

% Q = population growth that is used to "scale" resource extraction
% (variables other than fossil fuel exploitation are in per capita terms)
% popgr is the exogenous population growth rate

Q0 = 1;
popgr = 0.01;

% gamma is the coefficient of relative rsik aversion. If we wish
% to calibrate to a particular initial consumption level, we can
% allow gamma to vary.
gamma = 4.0;
gS_0 = 0.00015;  % Partial derivative of g with respect to S at t=0
% gamma = 4;
% gS_0 = 0.0001,

% Note that parms is a 1x14 vector

% We have three variables to choose inorder to hit three targets.


%%%%%%%%%%%%%%%  Regime 1: B>0, j>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

lambda0 = c0^(-gamma);
alpha1 = gS_0.*Res_0.^2;
alpha0 = g_0 - alpha1./Res_0;

parms = [delta A Sbar alpha0 alpha1 alpha2 alpha3 Gamma1 H_Tbar alpha psi beta gamma Q0 popgr Abar];

%======================================
% 1. Guess the time Tbar when the economy transitions to the analytical regime:
% Tbar = 321.97904698;
% k_Tbar = 1.417049897989985e+005;
% S_T0 = 1.640303327183687e+003; 
Tbar = 315.8;
k_Tbar = 141704.98998437249;
S_T0 = 1613; 
% Tbar = 500;
% k_Tbar = 108544.121458981;
% S_T0 = 1620.2025327183;



disp(' ');
disp(['Guessed values: Tbar = ' num2str(Tbar,15) ' k_Tbar = ' num2str(k_Tbar,20) ' S_T0 = ' num2str(S_T0,15)]);

%=======================================


% At Tbar, the marginal cost of backstop technology becomes Gamma2. Hence
% we have (Gamma1+H)^(-alpha)=Gamma2

H_Tbar = Gamma2^(-1/alpha)-Gamma1;

% The shadow value, eta, of H at T_bar must also equal zero since further changes
% in H have no effect on maximized utility beyond T_bar

eta_Tbar = 0;

%======================================

% Given Tbar and k_Tbar we can solve for Kbar, that is, the constant in the analytical solutions
% for k and lambda in the terminal regime.

Kbar = (k_Tbar*(beta*gamma-Abar*(gamma-1))*exp(Abar*Tbar/gamma)/gamma)^(-gamma);

% Having solved for Kbar, the value of lambda at the beginning of the analytical regime then
% is determined:

lambda_Tbar = Kbar*exp(Abar*Tbar);

%======================================================
% Now we have the initial conditions to solve the renewable regime backward to the
% time T0 when H=0, which represents the transition point from fossil fuels to renewables

z_Tbar=[k_Tbar,H_Tbar,lambda_Tbar,eta_Tbar];
tspan1 = [Tbar 0];

% events1 detects H = 0

options = odeset('RelTol',5e-14,'AbsTol',5e-14,'events',@events1);
[tr,zr,TE1,ZE1,IE1] = ode45(@renew,tspan1,z_Tbar,options);

% Extract components of z for graphing

kr=zr(:,1);
Hr=zr(:,2);
lambdar=zr(:,3);
etar=zr(:,4);

% Calculate p, c, j and inv (investment in k) and epsilon in the renewable regime for graphing

pr = (Gamma1+Hr).^(-alpha);
cr = lambdar.^(-1/gamma);
jr = A.*kr.*(etar.*(1-psi)./lambdar).^(1/psi);
invr = A.*kr.*(1-pr)-cr-jr;
epsrelr = pr - psi.*(1-psi).^(1./psi-1).*(lambdar.^(-1/psi)).*(etar.^(1./psi));

figure(1)

subplot(2,4,1), plot(tr,kr)
title('Capital k - renewable regime')
subplot(2,4,2), plot(tr,Hr)
title('Cumulative knowledge - renewable regime')
subplot(2,4,3), plot(tr,etar)
title('Shadow price of H - eta')
subplot(2,4,4), plot(tr,lambdar)
title('Shadow price lambda - renewable regime')
subplot(2,4,5), plot(tr,invr)
title('Investment - renewable regime')
subplot(2,4,6), plot(tr,jr)
title('Technology investment j')
subplot(2,4,7), plot(tr,pr)
title('Renewable energy production cost')
subplot(2,4,8), plot(tr,epsrelr)
title('Energy relative price eps/lambda')

% Extract initial values to use as terminal values in the fossil regime.

T0 = TE1;
k_T0 = ZE1(1);
lambda_T0 = ZE1(3);
eta_T0 = ZE1(4);

%%%%%%%%%%%%%%%  Regime 0: R>0, n>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Use the shadow values at T0 to specify starting values for sigma and nu

sigma_T0 = 0;
nu_T0 = 0;

%======================================


% Use equality of the shadow price of energy at T0 to solve for N0 given S0

eps_T0 = lambda_T0*Gamma1^(-alpha) - psi*((1-psi)^(1/psi-1))*...
    (eta_T0^(1/psi))*(lambda_T0^(1-1/psi));

g_T0 = eps_T0/lambda_T0;
N_T0 = alpha2/(Sbar-S_T0-alpha1/(g_T0-alpha0))-alpha3;

z_T0 = [k_T0,S_T0,N_T0,lambda_T0,sigma_T0,nu_T0];

tspan0 = [T0 0];
% %
% options = odeset('RelTol',1e-6,'AbsTol',1e-6,'events',@events0);
% [tf,zf,TE0,ZE0,IE0] = ode45(@foss,tspan0,z_T0,options);
% 
options = odeset('RelTol',5e-14,'AbsTol',5e-14);
[tf,zf] = ode45(@foss,tspan0,z_T0,options);

% Extract components of z for graphing

kf=zf(:,1);
Sf=zf(:,2);
Nf=zf(:,3);
lambdaf=zf(:,4);
sigmaf=zf(:,5);
nuf=zf(:,6);

% Calculate c, g and the derivatives of g for graphing

cf = lambdaf.^(-1/gamma);
g = alpha0 + alpha1./(Sbar-alpha2./(alpha3+Nf)-Sf);
gdS = alpha1*(alpha3+Nf).^2./((Sbar-Sf).*(alpha3+Nf)-alpha2).^2;
gdN = -alpha1.*alpha2./((Sbar-Sf).*(alpha3+Nf)-alpha2).^2;
gd2N = 2.*alpha1.*alpha2.*(Sbar-Sf)./((Sbar-Sf).*(alpha3+Nf)-alpha2).^3;
gdSN = -2*alpha1.*alpha2.*(alpha3+Nf)./((Sbar-Sf).*(alpha3+Nf)-alpha2).^3;


% Calculate n, inv and epsilon for graphing

Qf = Q0.*exp(popgr*tf);
nzero = (lambdaf > nuf);

% T1 = sum(nzero)+1;
% cond43 = (delta+g(T1).*A-gdN(T1).*A.*kf(T1)-A).*lambdaf(T1) - sigmaf(T1).*Qf(T1).*A;

npos = (gdN.*kf.*(delta+g.*A-A+sigmaf.*A.*Qf./lambdaf)+cf.*gdN-gdSN.*Qf.*A.*kf.^2-...
		popgr.*sigmaf.*Qf./lambdaf)./(gd2N.*kf-2.*gdN);

nf = npos.*(ones(size(tf))-nzero);
invf = A.*kf.*(1-g)-cf-nf;
epsrelf = g - sigmaf.*Qf./lambdaf;

% % Extract the values of S, N and lambda at TE0 so we can iterate to convergence
% 
% S_0 = ZE0(2);
% N_0 = ZE0(3);
% lambda_0 = ZE0(4);
% 
% % Compare calculated values to the target values
% 
% disp(['These should be zero: TE0 = ' num2str(TE0) ' S_0 = ' num2str(S_0) ' N_0 = ' num2str(N_0)]);
% disp(['Calculated lambda_0 = ' num2str(lambda_0) ' Desired value = ' num2str(lambda0)]);

% Extract the values of S, N and lambda at TE0 so we can iterate to convergence
k_0 = zf(length(tf),1);
S_0 = zf(length(tf),2);
N_0 = zf(length(tf),3);
lambda_0 = zf(length(tf),4);

k0_diff = k_0-k0;
% Compare calculated values to the target values
dX = norm([k0_diff S_0 N_0]);

disp(['These should be zero: k0_diff = ' num2str(k0_diff,10) ' S_0 = ' num2str(S_0,10) ' N_0 = ' num2str(N_0,10) ' dX = ' num2str(dX,8)]);

c0_calc = lambda_0^(-1/gamma);

disp(['Calculated c_0 = ' num2str(c0_calc) ' Desired value = ' num2str(c0)]);

n0_calc = nf(length(tf));

disp(['Calculated n_0 = ' num2str(n0_calc) ' Desired value = ' num2str(n0)]);

figure(2)

subplot(3,3,1), plot(tf,kf)
title('Capital k - fossil regime')
subplot(3,3,2), plot(tf,Sf)
title('Fuel use S - fossil regime')
subplot(3,3,3), plot(tf,Nf)
title('Mining technology capital N')
subplot(3,3,4), plot(tf,lambdaf)
title('Shadow price lambda - fossil regime')
subplot(3,3,5), plot(tf,sigmaf)
title('Fossil fuel user cost sigma')
subplot(3,3,6), plot(tf,nf)
title('Mining technology investment n')
subplot(3,3,7), plot(tf,invf)
title('Capital investment')
subplot(3,3,8), plot(tf,g)
title('Mining cost g')
subplot(3,3,9), plot(tf,epsrelf)
title('Energy relative price eps/lambda')

yr = A*kr;
yf = A*kf;
y = [yr;yf];
costf = g.* yf;
costr = pr.*yr;

figure(3)

subplot(2,3,1), plot([tr;tf],[kr;kf])
title('Capital k')
subplot(2,3,2), plot([tr;tf],[cr;cf])
title('Consumption c')
subplot(2,3,3), plot([tr;tf],[invr;invf])
title('Capital investment i')
subplot(2,3,4), plot([tr;tf],y)
title('Total Output Ak')
subplot(2,3,5), plot([tr;tf],[lambdar;lambdaf])
title('Shadow price lambda')
subplot(2,3,6), plot([tr;tf],[jr;nf])
title('Investment in Technology n&j')

figure(4)
subplot(2,3,1), plot([tr;tf],[cr;cf]./y)
title('Consumption Share of Output c/Ak')
subplot(2,3,2), plot([tr;tf],[epsrelr;epsrelf])
title('Energy Relative Price eps/lambda')
subplot(2,3,3), 
plot([tr;tf],[invr;invf]./y)
title('Capital Inv Share of Output i/Ak')
subplot(2,3,4), plot([tr;tf],[jr./yr;nf./yf])
title('Tech Inv Share of Output n&j/Ak')
subplot(2,3,5), 
plot([tr;tf],[(invr+jr);(invf+nf)]./y)
title('Total Inv Share of Output (i+n&j)/Ak')
subplot(2,3,6), plot([tr;tf],[costr;costf]./y)
title('Energy Mining or Production Cost Share of Output (gR+pB)/Ak')
